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ISOTHERMAL  COMPRESSIBILITY  OF  SPC/E  WATER 

Kazi  A.  Motakabbir  and  M.  Berkowitz 

Department  of  Chemistry,  University  of  North  Carohna  at  Chapel  Hill, 

Chapel  Hill,  North  Carolina,  27599 

Abstract 

Molecular  dynamics  computer  simulations  on  rigid  SPC/E  water  molecules  were 
performed.  The  goal  of  the  simulations  is  to  study  the  behavior  of  isothermal 
compressibility,  which  was  calculated  in  the  simulations  at  T=298,  273  and  248K  The 
calculated  isothermal  compressibilities  at  these  temperatmes  display  a  trend  contrary  to 
the  experimental  observations.  The  hydrogen  bonded  network  in  water  was  also 
investigated.  No  correlation  between  the  temperature  dependence  of  isothermal 
compressibility  and  the  munber  of  hydrogen  bonded  pentagons  was  observed. 
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I.  Introduction 


Among  many  unusual  features  displayed  by  liquid  water  the  anomalous  behavior  of 
thermodynamic  response  functions  is  one  of  the  most  interesting.  The  anomaly  in  the 
behavior  of  the  response  functions  such  as  constant  pressure  heat  capacity  (Cp),  thermal 
expansivity  coefficient  (a^)  and  isothermal  compressibility  («:^),  especially  at  temperatures 
below  0°C  is  believed  to  be  due  to  the  onset  of  some  sort  of  cooperative  phenomena, 
implying  the  possible  presence  of  anomalous  fluctuations  in  density  [1],  Since  the 
isothermal  compressibility  is  a  direct  measure  of  density  fluctuations  we  shall 
concentrate  our  attention  on  this  quantity.  While  for  most  of  the  liquids  k,^  decreases 
when  the  temperature  decreases,  for  liquid  water  k^,  displays  a  minimum  at  T  «  50°C. 
Below  this  temperature  rises,  especially  very  rapidly  in  the  supercooled  region.  The 
described  anomalous  behavior  of  the  compressibility  is  attributed  to  the  superposition  of 
"normal"  fluctuations  and  anomalous  fluctuations,  that  diverge  at  temperature  Tg  [2].  The 
nature  of  these  anomalous  fluctuations  and  their  divergence  at  Tg  is  not  clear. 

Several  theoretical  models  try  to  explain  these  anomalous  density  fluctuations. 
Stanley  and  Teixeira  [3]  ascribe  the  cooperative  nature  of  the  anomalous  fluctuations  to 
the  clustering  of  forufold  hydrogen  bonded  water  molecules.  The  patches  of  such  clusters 
are  characterized  by  local  order  and  reduced  local  density.  The  Stanley-Teixeira  model  is 
a  percolation  model  with  the  percolation  threshold  for  fourfolded  molecules  achieved  at  Tg. 
Another  model  is  due  to  Stillinger  [4],  who  proposed  the  existence  of  a  significant 
concentration  of  bulky  imstrained  polyhedra  in  water.  The  cooperativity  of  the  fluctuations 
is  again  due  to  the  geometry;  for  example  a  pair  of  dodecahedra  can  be  stabilized  by 
sharing  the  common  pentagonal  face.  The  model  of  Speedy  [5]  is  similar  to  Stillinger’s. 
While  Stilhnger  does  not  specify  the  particular  polygonal  structure,  Speedy  concentrates 
his  attention  on  pentagonal  hydrogen  bonded  rings.  These  rings,  according  to  Speedy,  have 
a  tendency  to  self  replicate,  thus  generating  low  density  patches.  In  Speedy’s  model  the 
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density  fluctuations  grow  proportionally  to  the  number  of  five-membered  hydrogen  bonded 
rings  until  the  limit  of  stability  is  approached  at  Tg. 

From  the  outset,  computer  simulations  helped  us  to  choose  among  the  competing 
theories,  proposed  to  explain  the  experiment.  Thus  the  first  molecular  dynamics 
simulation  of  water  by  Rahman  and  Stillinger  [6]  clearly  displayed  that  water  is  correctly 
described  by  a  distorted  hydrogen  bonding  network  or  "continuum"  model  and  not  by  its 
rival  the  mixture/interstitial  model.  One  may  also  hope  that  computer  simulations  will 
assist  us  in  determining  which  model  among  the  three  described  above  can  explmn  the 
nature  of  anomalous  fluctuations  in  density.  Indeed,  simulations  that  address  the  question 
of  the  nature  of  density  fluctuations  at  low  temperature  are  found  in  the  literature.  The 
results  are  controversial  at  least.  To  support  Stanley  and  Teixeira’s  theory  Geiger  and 
Stanley  performed  a  MD  simulation  on  ST2  water  [7],  in  which  they  observed  the  existence 
of  low-density  patches  of  water  molecules  [8].  Rapaport  performed  a  MD  simulation  of 
MCY  water  [9  ]  at  -22‘’C  and  did  not  observe  any  correlation  between  local  density  and  the 
number  of  hydrogen  bonds  associated  with  each  water  molecule  [10].  Speedy  and  Mezei 
in  their  work  aimed  to  support  Speedy’s  theory  reported  an  increase  in  five-member  rings 
concentration  with  the  temperature  decrease  [11].  Therefore  they  concluded  that  water 
anomalies  may  be  related  to  self-replicating  propeiisity  of  pentagons.  The  problem  with 
their  analysis  is  that  they  have  correlated  the  increase  in  pentagon  concentration  observed 
in  "computer"  water  with  the  increase  in  the  compressibility  in  "real"  water.  But  what  if 
the  compressibility  of  the  "computer"  water  (which  of  course  may  depend  on  the  model 
used  in  the  simulation)  is  not  increasing  when  temperature  is  decreased?  To  reach  defmite 
conclusions  on  the  correlation  between  the  polygon  concentration  increase  and 
compressibility  increase,  one  has  to  calculate  the  compressibility  of  the  "computer"  water. 
We  know  of  only  one  study  where  the  direct  calculation  of  compessibility  and  its 
temperature  dependence  was  performed.  This  was  done  by  Jorgensen  and  Madura  in  their 
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Monte  Carlo  simulation  performed  on  TIP4P  potential  [12],  The  reported  values  for  /c-j.  in 

[12]  are:  67  ±  13  x  10' '  Pa'  at  T  =  25°C  (the  experimental  value  is  45.2  x  10"  Pa'),  24 

±  3  X  10"  Pa  '  at  T  =  0°C  (50.9  x  10"  Pa  '  in  experiment)  and  23  ±  3  x  10"  Pa  '  at  T  =  - 

25°C  (70.9  X  10"  Pa  '  in  experiment).  One  can  see  that  according  to  ref.  [12],  for 

TIP4P  water  does  not  behave  in  the  right  way.  The  wrong  behavior  is  also  observed  for 

Cp.  The  problem  may  be  with  TIP4P  potential,  but  most  probably  it  is  in  the  length  of 

the  simulation.  The  simulation  was  performed  under  constant  pressure  and  temperature 

conditions  and  the  values  of  response  functions  Cp  and  x-j.  were  calculated  using  fluctuation 

formulas.  To  get  good  quantitative  numbers  using  fluctuation  formulas  one  needs  very  long 

simulations.  In  this  respect  calculations  of  Cp  or  are  similar  to  the  calculation  of 

dielectric  constant  e,  which  is  known  to  require  a  very  long  simulation  time  [13].  To  avoid 

very  long  simulations  one  can  use  the  thermodynamic  definition  of  the  response  function. 

In  this  case,  to  find  from  the  computer  experiment  one  has  to  determine  that  portion 

of  the  isotherm  where  one  is  interested  in  its  slope.  That  can  be  obtained  from  2  or  3 

relatively  short  simulations.  In  what  follows  we  describe  such  simulations  and  present  the 

values  of  isothermal  compressibility  for  SPC/E  water  [14]  at  three  different  temperatures. 

We  also  discuss  the  connection  between  the  compressibility  behavior  and  the  concentration 
/ 

of  five  membered  hydrogen-bonded  rings  in  water  structure. 

II.  Molecular  Djniamics  Simulations 

Six  molecular  dynamics  simulations  in  a  canonical  N,V,T  ensemble  and  one 
simulation  in  N,P,T  ensemble  were  performed.  For  N,V,T  ensemble  two  simulations  were 
performed  at  T=248  K,  two  simulations  at  T=273  K  and  two  simulations  at  T=298  K  For 
the  canonical  ensemble  simulations  at  every  temperature  one  simulation  was  performed  at 
the  density  roughly  corresponding  to  1  g/cc  and  at  pressure  aroxmd  1  atm.,  and  the  second 
simulation  at  the  elevated  density  and  elevated  pressure.  Every  simulation  was  of  60  ps 
duration  time,  20  ps  were  used  for  equilibration  and  the  remaining  40  ps  were  used  for 
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analysis.  The  40  ps  were  divided  into  two  blocks  of  20  ps  each  and  the  average  pressure 
over  this  block  of  time  was  calculated.  Only  one  N,P,T  simulation  was  performed  at  T  =  248 
K.  This  simulation  was  done  for  longer  duration,  400  ps,  so  that  we  can  use  the 
fluctuation  formula  for  compressibility  calculation. 

The  maintenance  of  a  constant  temperature  and/or  pressure  in  the  system  was 
achieved  by  an  application  of  a  Ferrario-Ryckaert  algorithm  [15],  based  on  the  extended 
variable  method  of  Andersen  [16]  and  Nose  [17].  To  find  the  values  of  the  inertial 
parameters  Wg  and  Wq  that  drive  the  fluctuations  in  temperature  and  voh  ne,  several  trial 
simulations  were  performed.  We  finally  settled  for  values  Wg  =  100  (kJ/mol)ps^  and  Wq 
=250  (kJ/mol/nm®)ps^,  since  these  were  the  values  that  produced  the  most  stable 
molecular  d3mamics.  The  potential  energy  of  interaction  between  two  water  molecules  was 
multiplied  by  the  switching  fimction  proposed  by  Steinhauser  [18],  which  smoothly  reduces 
the  energy  from  its  value  at  Rr= 0.8075  nm  to  zero  at  Rc  =  0.85  nm.,  the  long-range  forces 
contribution  was  taken  into  account  by  usage  of  the  reaction  field  method.  The  dynamics 
was  carried  out  on  216  molecules  with  the  use  of  the  Verlet  algorithm  [19]  and  the  SHAKE 
procedure  [20].  The  time  step  was  2.5  fs  and  every  fourth  configuration  was  saved.  The 
water  model  used  in  the  simulations  was  SPC/E  model.  It  was  recently  shown  that  this 
model  in  a  nice  way  reproduces  many  of  the  features  of  bulk  water  [14,21]. 

III.  Results  and  Discussion 

a.  Compressibility  Calculation 

The  canonical  N,V,T  simulations  were  very  similar  to  the  ones  recently  performed 
in  the  study  of  pressme  dependent  properties  of  TIP4P  water  [22].  As  in  the  simulation 
with  T1P4P  water  the  present  simulation  displays  the  weakening  of  the  tetrahedral  order 
taking  place  in  water  with  the  density  increase.  This  is  revealed  by  a  diminished  intensity 
of  the  first  peak  and  by  gradual  disappearance  of  the  second  peak  with  the  density 
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increase,  illustrated  in  figure  1  for  T  =  248  K. 

To  calculate  the  compressibility  K-p  we  used  the  finite  difference  approximation,  i.e. 


The  linear  approximation  in  (1)  is  expected  to  be  good  for  water  in  the  region  0-2000 
atmospheres,  since  the  isotherm  for  real  water  is  nearly  linear  in  this  region.  In  equation 
(1)  Pi  is  the  density  corresponding  to  pressure  P,  and  pj  is  the  density  corresponding  to 
P2.  These  values  of  p  and  P  were  obtained  from  molecular  dynamics  runs  af  three 
different  temperatures.  To  enrich  the  statistics  every  run  was  divided  into  two  blocks  and 
for  every  block  the  values  of  P  were  calculated.  This  helps  to  explain  the  notation  used 
in  Table  1.  Thus  P,,  means  pressure  P,  corresponding  to  density  pi,  P,2  is  the  pressure 
corresponding  to  density  p,  in  the  second  half  of  the  run,  P21  is  the  pressure  corresponding 
to  density  p2  from  the  first  half  of  the  rim,  etc.  Such  scheme  permits  us  to  calculate  four 
possible  values  for  the  isothermal  compressibility  at  each  temperature.  The  values  of 
these  four  quantities  are  given  in  Table  1  together  with  the  average  value. 

We  have  also  used  the  N,P,T  ensemble  to  calculate  the  compressibility  of  water. 
Since  this  is  done  through  the  fluctuation  formula,  =  <AV2>/Vk3T,  it  demands  long 
simulation  times.  We  watched  the  time  dependence  of  the  cumulative  value  of  and 
when  it  reached  a  plateau  value  at  «35.5x1(1"  Pa  *  after  around  400  ps  the  simulation 
was  terminated.  The  cumulative  behaviour  of  Ktj,  is  shown  in  Fig.  2.  To  obtain  an  idea 
of  the  error  in  x-j,  present  in  this  simulation  we  divided  400  ps  into  three  blocks  200  ps 
each  (0-200ps,  100-300ps  and  200-400ps).  The  values  of  of  /c.j,  for  each  block  and  the 
average  value  are  given  in  Table  1.  As  one  can  see  the  values  of  K.p  obtained  from  N,P,T 
and  N,V,T  simulations  show  rather  good  agreement  between  themselves.  The  calculated 
values  for  the  isothermal  compressibility  display  a  slight  minimum  around  T«273K,  similar 
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to  the  one  found  experimentally  at  T«323K.  This  behavior  is  contrary  to  experimental 
observation  that  isothermal  compressibility  monotonically  increases  and  nearly  doubles  its 
value  when  the  temperature  is  lowered  from  T  =  298K  to  T  =  248K  It  is  therefore 
possible  that  the  thermodynamic  properties  of  a  "computer"  water  are  shifted  down  in 
temperature.  Such  downshift  was  for  example  foimd  for  TIP4P  water,  which  is  estimated 
to  have  a  melting  temperature  somewhere  around  240K  [24].  We  do  not  know  what  is  the 
melting  point  of  SPC/E  water,  but  it  is  possible  that  it  is  also  shifted  down  in  temperature. 

To  see  if  there  is  any  correlation  between  the  observed  behavior  of  x-j,  and  the 
concentration  of  pentagons  in  water  we  performed  a  topological  analysis  of  SPC/E  water 
structure,  which  we  describe  below. 

b.  Hydrogen  Bonds  and  Their  Network 

To  investigate  tmderl3dng  geometrical  structures  in  further  detail,  we  have 
enumerated  distribution  of  hydrogen  bonds  and  non-short-circuited  polygons  formed  by  the 
hydrogen  bonds.  Our  calculation  is  based  on  a  set  of  100  equidistant  configurations  taken 
from  the  final  20ps  run  at  lower  density  (/>,)  at  each  temperature.  In  our  simulation  a  pair 
of  molecules  is  considered  hydrogen  bonded  if  they  are  within  0.32  nm  away  (roughly  first 
minimum  of  the  oxygen-03Q^gen  radial  distribution)  and  their  interaction  energy  falls  below 
the  cutoff,  Vhb  =  -2.5  Kcal/Mole.  The  results  for  the  hydrogen  bond  distribution  are 
displayed  in  Table  2.  The  distribution  at  T=298  K  agrees  quite  well  with  previous  findings 
[23].  As  expected,  the  concentration  of  four-bonded  molecules  considerably  increases  as 
temperature  is  decreeised. 

To  calculate  the  ring  distribution  we  have  followed  the  definition  of  non-short- 
circuited  polygons  as  described  by  Rahman  and  Stillinger  [25],  viz,  no  two  vertices  of 
polygon  are  connected  by  shorter  number  of  bonds  than  the  miniminn  number  of  bonds 
already  existing  within  the  polygon.  In  conformity  with  their  concern  about  the  erronious 
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polygon  closures  due  to  the  small  number  of  molecules  (216)  used  in  the  simulation,  we 
have  used  eight  replicas  of  the  MD  cells  to  construct  a  2x2x2  larger  cube  containing  1728 
distinct  molecules  under  periodic  boundary  conditions  for  the  purpose  of  the  polygon 
enumeration. 

We  have  used  the  following  polygon  enumeration  algorithm  which  is  similar  to  the 
triplet  counting  algorithm  developed  by  Belch  and  Rice  [26].  In  this  scheme  for  a  given 
configuration  we  first  generate  a  list  of  bonded  neighbors  for  each  molecule.  A  molecule 
is  considered  for  polygon  search  only  if  it  has  at  least  two  bonded  neighbors.  For  &  given 
molecule  designated  as  ’B’  we  arbitrarily  designate  its  two  neighbors  as  ’A’  and  ’C’.  We 
start  the  search  by  looking  into  the  neighbors  of  ’C’.  If  ’A’  is  a  neighbor  of  ’C’  then  we 
conclude  ABC  is  a  triangle.  If  triangles  are  not  found  then  subsequent  neighbors  are 
searched  to  find  polygons  of  larger  sizes.  For  each  triplet  of  molecule  only  the  shortest 
non-short-circuited  polygon  is  listed.  After  exhaustively  searching  through  all  the 
molecules  using  all  distinct  triplets  and  enumerating  all  the  polygons  of  different  sizes  a 
final  search  is  performed  to  identify  duplicates  as  well  as  short  circuiting  of  larger  polygons 
by  smaller  ones. 

Using  the  above  scheme,  we  have  exhaustively  identified  all  possible  non-short- 
circuited  polygons  of  sizes  3  through  11.  We  did  not  look  for  larger  size  polygons  since  it 
would  consume  prohibitive  computer  time  and  also  not  very  pertaining  to  our  present 
study.  In  Table  3,  we  show  the  distribution  of  hydrogen  bonded  polygons  at  three 
temperatures  we  have  studied.  Although  pentagons  increases  in  number  with  temperature 
decrease  the  predominance  of  hexagons  becomes  more  pronounced  at  lower  temperature 
indicating  a  competition  between  nucleation  process  and  retention  of  hquid  structure. 
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IV  Conclusions 

Our  study  of  hydrogen  bonded  network  in  SPC/E  water  demonstrates  that 
temperature  decrease  does  not  produce  a  structure  where  five-membered  rings  are 
dominant.  As  a  matter  of  fact,  six  membered  rings  remain  dominant  and  the  rate  of  then- 
increase  is  parallel  to  the  rate  of  the  increase  in  number  of  five-membered  rings.  Similar 
behavior  was  recently  observed  in  the  simulation  performed  on  ST2  water  by  Mausbach, 
Schnitker  and  Geiger  [27].  We  agree  with  their  conclusion  that  the  observed  structural 
changes  due  to  the  temperature  decrease  can  be  attributed  to  the  approach -of  the 
amorphous  structure.  In  this  structure  hydrogen  bonded  rings  of  all  sizes  are  expected. 
Moreover,  a  lack  of  correlation  between  the  behavior  of  the  calculated  isothermal 
compressibility  and  the  nrunber  of  five-membered  rings  indicates  that  Speedy’s  hypothesis 
may  not  be  right.  We  also  would  like  to  point  out,  that  while  an  effective  pair-potential 
used  in  the  simulations  with  SPC/E  water  is  rather  successful  in  reproducing  many 
properties  of  liquid  water,  it  fails  to  reproduce  the  anomalous  behavior  of  isothermal 
compressibility  in  temperature  range  298K-248K  It  would  be  interesting  to  see  if  polar- 
polarizable  models  of  water,  which  are  now  under  investigations  [28,29],  demonstrate  the 
anomalous  behavior  of  thermodynamic  response  functions  in  this  temperature  range.  Note 
that  recently  designed  random  network  models  also  are  not  able  to  reproduce  the 
anomalous  behavior  of  the  isothermal  compressibility  [30,31].  We  hope  that  future 
computer  studies  of  low  temperature  behavior  of  water  will  help  to  resolve  the  natme  of 
the  anomalous  fluctuations  responsible  for  anomalous  behavior  of  the  response  functions. 
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Legends  for  Tables  ; 


1.  Isothermal  compressibility  of  SPC/E  water.  Density  (p  in  gm/cc),  pressure  (P 
in  atm)  and  compressibility  (/c  in  10"  Pa')  calculated  from  20  ps  blocks  are  shown  for 
NVT  ensemble.  Compressibility  using  NPT  ensemble  at  248  K  is  also  shown. 


2.  Distribution  of  hydrogen  bonds  per  molecule  in  SPC/E  water  at  three  different 
temperatures.  The  average  number  of  hydrogen  bonds  per  molecule  is  given  in  the  last 
row. 


3.  Distribution  of  hydrogen  bonded  polygons  in  SPC/E  water  and  their  temperature 
dependence.  The  actual  number  of  rings  alongwith  their  uncertainties  are  shown. 
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Legends  for  Figures  ; 


1.  The  effect  of  pressure  on  oxygen-oxygen  radial  distribution  at  248  K  in  NVT 
ensemble.  The  solid  line  is  at  lower  density  (1.00722  g/cc)  and  the  dotted  line  is  at  higher 
density  (1.0822  g/cc). 

2.  Cumulative  plot  of  isothermal  compressibility  versus  time  at  248  K  in  NPT 
ensemble. 
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TABLE  1 


T  248K 


273K  298K 


1.0072  1.0055  1,0033 

25.39  -280.42  278.08 

-68.58  -30.27  207.94 


P2  1.0822  1.0755  1.0550 


1779.97 

1762.16 

1450.79 

P22 

1790.14 

1698.72 

1385.60 

/Cl, 

39.58 

31.86 

41.43 

/C12 

39.35 

32.88 

43.87 

/C21 

37.57 

36.31 

39.05 

/C22 

37.36 

37.63 

41.25 

/Ct(NVT) 

38.47  ± 

1.16 

34.67  ±  2.74 

41.41  +  1.95 

/Ct(NPT) 

34.07  ± 

5.21 
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TABLE  2 


248K 

273K 

298K 

0 

0.00005 

0.0001 

0.00097 

1 

0.0038 

0.0088 

0.0153 

2 

0.0546 

0.0922 

0.1198 

3 

0.2918 

0.3486 

0.3859 

4 

0.6091 

0.5072 

0.4362 

5 

0.0403 

0.0425 

0.0411 

6 

0.0003 

0.0006 

0.0007 

<nb> 

3.63 

3.48 

3.37 
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TABLE  3 


Nj(T) 


j 

248K 

273K 

298K 

3 

2.1  ±  2.22 

2.6  ±  2.22 

3.3  ±  2.38 

4 

19.7  ±  4.78 

20.7  ±  5.03 

18.4  ±  5.31 

5 

80.2  ±12.62 

62.6  ±10.40 

54.9  ±11.09 

6 

106.7  ±14.47 

79.2  ±13.96 

62.5  ±11.92 

7 

87.6  ±12.26 

77.4  ±15.23 

61.4  ±12.15 

8 

52.6  ±  9.45 

49.1  ±  9.78 

47.0  ±10.25 

9 

26.1  ±  7.70 

31.8  ±  7.35 

33.5  ±  7.35 

10 

8.1  ±  4.67 

15.2  ±  5.62 

20.1  ±  5.51 

11 

2.0  ±  2.21 

6.1  ±  3.60 

9.5  ±  4.62 
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